mtest = matrix(3,3,i,j,i-j);

\\mu is global variable which must have been created before.
grams(m) = {
for(i=2,length(m),
  for(j=1,i-1,  mu[i,j] = (m[,j]~*m[,i]) / norml2(m[,j]);
    m[,i] -= m[,j] * mu[i,j])
);
return(m);
}


islll(x,lc) ={ \\ test if it is lll reduced, return 1 if it is the case.
my(tmp,n,l,r);
n = length(x);
mu = matrix(n,n);
x = grams(x); 
for(i=2,n, 
  if(l = norml2(x[,i] + mu[i,i-1]*x[,i-1]) < r = lc * norml2(x[,i-1]),
    print("1e\t",l,"\t",r,"\t",i,mu[i,i-1]); return(0) );
);

for(i=1,n,
  for(j=2,i-1, if( abs(mu[i,j]) > 0.5, print("2e"); return(0) ))
);
return(1);
}


Knapsack_generate(size,M) = { \\ generate Knapsack instances
my(a,x,tc,maxx);
x = vector(size);
a = vector(size);
t = 0;
for (i=1,size,
	x[i]=random(2); a[i]=random(M); t += x[i]*a[i]);
\\ get the density
maxx = vecsort(a,,4)[1];
d = size/log(maxx)*log(2);

tc = sum(i=1,size,a[i]) -t;
print("x = ", x); print("t = ", t); print("d = ", d); print("tc = ",tc); 
return(a);
}


Solve_Knapsack(a, t, num = 1) = {
\\ we work with the lines of the matrix, and do transpositions for feeding qflll
my(d,maxxn,m,l,j,lamb,tc,x);
n = length(a);
tc = sum(i=1,n,a[i]) -t; print("tc = ", tc);
\\ get the lattice from a and t
m = matid(n+1);
for(i=1,n, m[i,n+1] = a[i]);
m[n+1,n+1] = -t;

\\ give a lattice for fplll 
default(logfile, "matfplll" num);
default(log,1);
print("[");
for(i=1,n+1, print1("[");
  for(j=1,n+1,print1(m[i,j], " ")); print("]");
);
print("]");
default(log, 0);

l = m~ * qflll(m~,1);
print("Norme de l : ", norml2(l)); 
maxx = vecsort(a,,4)[1];
d = n/log(maxx)*log(2);
print("Density : ", d);

for(i=1,n+1, x = l~[i,];
  j = 1;
  while(x[j] == 0 && j <= n, j++);
  lamb = x[j];
  while((x[j] == 0 || x[j] == lamb) && j <= n, j++);
  if(j == n+1 && x[n+1] ==  0, x = x/lamb; 
    if(sum(i=1,n,x[i]*a[i]) == t, print("norme de x : ", norml2(l)); return(x))
  );
);

\\ 2nd lll, with complements
m[n+1,n+1] = -tc;
l = m~ * qflll(m~,1);
 
for(i=1,n+1, x = l~[i,];
  j = 1;
  while(x[j] == 0 && j <= n, j++);
  lamb = x[j];
  while((x[j] == 0 || x[j] == lamb) && j <= n, j++);
  if(j == n+1 , x = x/lamb; 
    if(sum(i=1,n, x[i]*a[i]) == tc, print("Complementaire! Norme de x : ", norml2(x)); 
      return(x)
    );
  );
);
\\ if failed, return the 2nd lll matrix
print("FAILED!!!  FAILED!!!  FAILED!!!  FAILED!!!  FAILED!!!  FAILED!!!");
return(norml2(l~));
}




